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1. Motivation 

The long-distance nature of QED poses a substantial challenge when including QED interac¬ 
tions in lattice QCD simulations. While it is possible to account analytically for the large power-law 
corrections that are introduced by the QED interaction in finite-volume QCD-i-QED simulations, 
see, e.g., [1, 2], such a setup requires the use of lattice ensembles that are much larger than nec¬ 
essary for a pure QCD simulation. An alternative procedure that explicitly couples regular QCD 
simulations in a finite volume (QCDy) to QED in infinite volume (QEDoo) is presented here. 

In the following we first introduce a specific prescription fo puf valence fermions and photons 
in infinite volume and fhen demonsfrafe fwo versions of fhe mefhod in a numerical applicafion 
to fhe leading-order hadronic vacuum polarizafion (HVP). We presenf resulfs for elecfro-magnefic 
mass splittings and give an ouflook fo general QCDy -i- QEDoc simulations. 

2. The setup 

Eet us consider a box of QCD gauge fields periodic in a volume V. We furthermore imagine 
that we repeat the QCD gauge field in each direction an infinite number of times and couple an 
infinite-volume valence fermion 'E to the repeated gauge background as illustrated in Eig. 1. The 
infinite-volume fermion is in turn coupled to infinite-volume photons via A^ {x)V^ (x), where 

is a vector cuiTent of fields *E. One may imagine that such a setup with photons living in infinite 
volume has largely suppressed QED finite-volume errors since there are no mirror charges and the 
finite-volume sum over photon momenta is explicitly replaced by the infinite-volume integral. 



Figure 1: Valence fermions 'E living on a repeated gluon background 1/^ with periodicity Ly L 2 and vectors 
Li = (Li,0),L2 = (0,T2). 

If restricted to a finite number of QCD gauge field copies, this setup is identical to solving 
the Dirac operator on the repeated QCD gauge background but with QED gauge fields living in 
the larger volume. In the following we present a prescription to generate the infinite-volume (or 
large-volume) setup stochastically. 
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Let 0® be a fermion field defined in finife volume wifh fwisfed boundary condifions 


A 


where jj. = {x,y,z,t}, Lf^ is fhe size of fhe QCD box in direcfion /i, ft is fhe unif vector in direction 
jj., and 6^ is the twist angle in direction /i. In analogy to Bloch’s theorem, the symmetry under 
translations by L^fi then yields 

pin ^4 a . , 

(<^’r ) > (^.l) 


where the (•) denotes the fermionic contraction in a fixed background gauge field U^{x), L = 
and x,y,n,m,6 are four-vectors. For a finife number of QCD gauge field copies in 
direcfion n, fhe infegral ^ is replaced by fhe sum (1 ^o. 2 ;r/iv^., 2 ;r(iv^- 1 ■ For = 

2 fhis reduces fo fhe well-known periodic-plus-anfiperiodic frick. The fwisf-average of Eq. (2.1) 
may be performed stochastically and can be used fo generafe fhe infinife-volume photon momenfum 
infegrals, see Sec. 4. 

The reduction of finife-volume errors fhrough averaging of solutions wifh differenf fwisfs is 
well-known in fhe sfudy of mefallic sysfems [3] buf was also recenfly investigated in the study of 
two-baryon systems [4] . 

For a general observable, we propose the following prescription: 


1. Use the infinite-volume fields 'T for all operafors and sources, 

2. perform fhe Wick confracfions, and finally 

3. use Eq. (2.1) fo wrife fhe resulf in ferms of infegrals over fwisfs involving only Dirac inver¬ 
sions of fhe finife-volume fheory. 


The posifions of fhe vertices in fhe resulting expression can fhen be summed over fhe infinife 
volume. For quark-connecfed diagrams such a sum can be performed naively as described below. 
For quark-disconnecfed confribufions one needs fo resfricf fhe position of separafe quark loops fo 
fhe same QCD box.^ Such a resfricfion can be readily implemenfed in fhe prescripfion discussed 
below. 

The sefup described in fhis secfion works well wifh multi-source mefhods such as AMA [5]. In 
fhe remainder of fhis manuscripf fhe numerical dafa uses an AMA sefup wifh only one fwist vecfor 
per source (such fhaf fhe number of fwisf vecfors fhaf are averaged coincides wifh fhe number of 
sources). Nofe fhaf one may be able fo re-use zero-fwisf eigenvecfors by employing solvers such 
as [6] fhaf use blocked eigenvecfors. 


3. The muon hadronic vacuum polarization 

The leading-order hadronic vacuum polarizafion diagram, shown in Fig. 2, is convenienfly 
evaluafed following [7] which expresses fhe full diagram as 

a^= d{q'^)F{q^)n{q^) (3.1) 

■Jo 

*We would like to thank Luchang Jin for bringing this point to our attention. 
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with function F{q^) defined in [7] capturing the photon and muon components of the diagram, 
n{q^) = _n('0), n^v(^) = {q^qv -q^ggv)'n{q^). Furthermore n^v(^) = 

and Ylfiy{x) = (V^(;c)Fv(0)) with vector current In the following numerical examples, we use 
a conserved vector current at position v and a local vector current at the origin. 



Figure 2: The leading-order hadronic vacuum polarization diagram. The lower fermion line is the muon. 

Equation (3.1) follows the prescription that we outlined above. We couple the QCD fermion 
loop to an infinite-volume photon whose propagators are included in the function F{q^). By com¬ 
puting n^v defined in terms of the infinite-volume fields 'F, we have access to a continuum of 
photon momenta q. The most significant simplification of this specific example is that we have re¬ 
duced the four-dimensional photon momentum integral analytically to a one-dimensional integral 
over q^. In this way, we can test our prescription in the case of twisting in just a single direction. 
In a general QCD-i-QED problem, we will evaluate four-dimensional photon-momentum integrals 
and appropriately included non-zero twist angles in all four directions. 

In the remainder of this section we first explore a naive implementation of the twist-averaging 
idea to evaluate a^. Eet C^v(0 = L.^n^vC-^o = C{t) = for /r = x,y,z. Then we can 

write down an estimator for tl{q^) that configuration-by-configuration satisfies = 0) = 0 and 
fhereby avoids a statistical noise problem for small momenta. Starting from = 0)) = 0 

and (lmn^v(?^)) = 0 we arrive at 

(fitf» = (^Re ReC„(r)) (3.2) 

which is a minimal modification of Eq. (81) of [8]. This expression includes an explicit double¬ 
subtraction of = 0) and Il{q^ = 0). 

More specifically, we compute 


X 



0 


averaging over independent temporal twist angles di and 02 for the two propagator lines. All spatial 
twists are zero such that our setup corresponds to an infinite time direction. We inject momentum 
only in time direction and cut the sum over time-slices at sufficiently large distance such that the 
contribution of the remainder can be ignored within the statistical precision. In Sec. 4, we present 
a method that does not require such a cut. 
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Figure 3 shows a numerical comparison of results using RBC/UKQCD’s 16^ and 24^ en¬ 
sembles that only differ in their respective physical volume. This allows for a test of remnant 
finite-volume errors introduced by the sea sector. The integral over is performed using both the 
Trapezoidal and Simpson’s rule, choosing a step-size such that their difference is below 1 /100 of 
the statistical error. We find that the relative error of is consistent with the almost -independent 
integrand uncertainty shown in Fig. 4. 
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Figure 3: Strange quark contribution on RBC/UKQCD’s 16^ and 24^ ensembles with « 1.73 GeV, 
mi = 0.01 (ifijc « 422 MeV), and = 0.04. The 16^ (24^) measurement was performed using 2 (1) exact 
and 32 (16) inexact sources on 60 (78) configurations. The left figure shows the times-lice t contribution to 
n(k = nip), where nip is the physical muon mass. The right figure shows the integrand of ap as a function 
of q^. Note that the stochastic noise is well-behaved for small momenta. 


4. Stochastic integration of photon momenta 


In this section we demonstrate an efficient method to use the twist-averaging procedure to 
sample over the photon momenta. A wise choice of sampling weight, i.e., the use of importance¬ 
sampling techniques can yield a substantial benefit. The following discussion is explicitly given 
in one dimension but all expressions and methods translate in a straightforward way to the more 
general four-dimensional case. 

We continue the discussion of the HVP diagram to illustrate the method. The full diagram 
with lattice regulator can be written as 



/ TT plK 

dkT>^^{k) / 

-n Jo 


^^2 y ^ik{nL+x)+in(ei 

In In ;cg{0,...,L-l},n6N 


-^)C^v(u0i,02) 


(4.1) 


for an appropriately defined F^'^. Poisson’s summation formula yields 

r dkr^^^ik) {62-ei)/L)Cpy{x,61,62) ( 4 . 2 ) 

J-^ Jo 271 271 L-I} 

with 5{k) = + By writing J^^dkg{k) = {Zn=o,...,L-i So""d^kg{J-T^n/L + 

6k/L) and performing the integral over 6k, we obtain a method to stochastically sample the photon 
momenta through the random choice of twist angles. 
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In Fig. 4 we show the resulting error for this method using the same configurations and sources 
as used in Fig. 3. We now, however, use importance sampling for the random twist angles such that 
they follow the probability distribution induced by of Eq. (3.1). We observe a slight reduction 
of statistical eiTor, which may demonstrate the benefit of the importance sampling strategy. 
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Figure 4: Relative errors of integrated up to maximal momentum (left) as well as the integrand of a^ 
for (right) for the methods described in this manuscript. 


The sum over L-translations in Eq. (4.1) generates conservation of the small-momentum com¬ 
ponent 6 /( = 62 — Oi- This leads to an intuitive pictorial prescription in momentum space 


— Zi<jpBC={0,±27r/i.,...} 

<7 pbc + ^qi L 



k + qpBC + ^q/L 


with average over twist-angle 6 g and sum over finite-volume periodic momenta ^pbc- 


5 . QED mass splittings 


A potentially interesting application of the ideas explained above is the computation of QED 
mass splittings. Eor a general meson two point function, one needs to evaluate (ignoring quark- 
disconnected diagrams for now) 

+ {JJ =Co{t) + aCi{t) 

for which the mass-shifts could be obtained from effective mass plots following 


meff(t) = meffp(t) -f (t) , 




Ci(t) Ci(f + 1) 
Co{t) Co(f + l)' 


In this general case a four-dimensional photon momentum integral and hence a four-dimensional 
twist-averaging procedure is necessary. Using an appropriate importance-sampling of twist angles, 
the methods outlined above should be applicable in a straightforward manner. The photon propa¬ 
gator could, e.g., be implemented stochastically in a sequential source setup. 

Eor a first test that allows us to re-use the HVP measurements, we use current algebra and soft 
pion theorems [9] and compute 
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Results of a quick numerical test are shown in Fig. 5. 
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Figure 5: V-A contribution for the strange quark loop. At large there is a gap due to the heavy mass 
such that the integral does not converge. The stochastic integration method of the photon momentum with 
{aqY < 1 cut yields an « 4% statistical error. 
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6. Conclusion 

We have proposed a prescription to put valence fermions and photons in infinite volume. We 
expect a substantial reduction of finite-volume errors in lattice QCD-i-QED simulations in this setup 
such that regular lattice QCD boxes could be used for combined QCD-i-QED simulations. We have 
discussed both a naive implementation of the idea and a refined sfrafegy fhaf samples over phofon 
momenfa sfochasfically. In fhe laffer mefhod an importance-sampling sfrafegy seems promising. 

In fhe fulure, fhe mefhods presenfed here will be fesfed for realistic pion masses in fhe confexf 
of QED mass spliffings, QED corrections fo decay consfanfs, as well as finile-volume correcfions 
fo fhe hadronic confribufions fo (g — 2)^, in particular fhe lighf-by-lighf confribufion. 
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